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FOREWORD 


This  work  reports  part  of  two  man-years  of  effort  supported  by  the 
Advanced  Lightweight  Torpedo  Project.  The  object  is  the  direct  prediction 
of  the  onset  and  details  of  submarine  pressure  hull  rupture  caused  by 
underwater  explosive  attack.  An  earlier  report  developed  a  material  model 
describing  plasticity  and  rupture  under  extremely  rapid  loading  conditions. 
This  work  gives  its  finite  element  implementation,  with  numerical  results 
to  be  communicated  in  a  subsequent  report. 
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INTRODUCTION 


In  ductile  structural  metals  such  as  aluminum,  the  response  to  rapidly 
applied  high  loads  can  involve  two  basic  mechanisms  [1,2]:  (a)  flow, 

understood  physically  as  slip  within  grains,  and  (b)  damage,  comprising  the 
nucleation  of  microvoids  at  grain  interfaces,  their  subsequent  growth,  and 
their  eventual,  usually  abrupt,  coalescence  into  cracks.  Reference  1 
introduces  a  constitutive  model  extending  viscoplasticity  to  accommodate 
both  flow  and  damage.  The  present  work  concerns  the  finite  element 
implementation  of  the  model.  Numerical  results  will  be  presented  in  a  later 
work. 


CONSTITUTIVE  MODEL 


We  briefly  state  the  constitutive  relations  given  in  Reference  1  under 
restriction  to  small  strains  and  temperature  independent  deformations. 

In  obvious  notation,  the  strain  is  decomposed  into  elastic,  flow  and 
damage  parts  according  to 


€ 


ij 


and  the  elastic  strain  is  given  by  Hooke's  law  as 


(1) 


3ij  ■ 2a  4 + '  4 


(2) 


where  u  and  \  are  the  Lame'  coefficients  and  5,.  is  the  Kronecker  tensor. 

1  J 

Let  e^.  and  e  be  the  deviatoric  (shear)  and  isotropic  (dilatational ) 
parts  of  the  strain  tensor,  and  s.,  and  s  correspondingly  for  In 

1  J  1  J 

Reference  1,  constitutive  relations  embodying  associated  flow  rules  were 


developed  as 

T.  Nicholson,  D.  W.,  "Constitutive  Model  for  Rapidly  Damaged  Structural 
Materials,"  accepted  for  publication  in  Acta  Mechanica 
2.  Barbee,  T.  W.,  et  al ,  "Oynamic  Fracture  Criteria  for  Ductile  and  Brittle 
Metals,"  J.  Materials,  1972 
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f  f  3Ff 

eij  '  "f  '  *f  (Ff  k  ^  >  9S .  . 

*  J 

(3a) 

-  0 

(3b) 

ed  -  nd  <  <d  (Fd  -  kd)  > 

(3c) 

4d.  -  0. 

(3d) 

Here  and  nd  are  material  constants,  kF  and  kd  are  parameters  representing 

dependence  on  the  history  of  flow  and  damage,  4>d>  Ff  and 

functions  and  the  symbols  <  •  >  are  defined  by 

F^  are  material 

o  »  r  i  o 

<  *  (r)  >  = 

r  ,  r  >  o  . 

The  material  function  F,  depends  on  e?.,  kF  and  s.  .,  while  F. 

1  ‘  J  1  J  u 

kd  and  s.  Finally,  the  history  parameters  are  governed  by 

depends  on  ed. 

kf  ■  "W*  ^  *pq>  <3 

(3e) 

kd  -  hd  (ed,  kd,  s)  ed 

(3f ) 

FINITE  ELEMENT  FORMULATION 


Equation  of  Equilibrium  for  an  Element 


Suppose  that  high  loads  are  raDidly  applied  to  a  body  governed  by 
Equation  3a-f.  The  body  can  be  represented  as  a  collection  of  finite  elements 
connected  to  each  other  at  nodes  [3].  We  consider  equilibrium  of  a  given 
element. 


In  accordance  with  the  usual  practice  in  finite  element  analysis,  we 


hereafter  use  vector  notation.  So  is  replaced  by  r,.  .  by  c_,  etc. 

T.  Zienkiewicz,  0.  C. ,  The  Finite  Element  Method,  Third  Edition,  McGraw-Hill 
Book  Co.,  New  York,  1977 
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Let  the  vector  r  denote  the  position  of  a  given  interior  point  of  the 
element  under  study.  The  time  displacement  vector  u_(r)  is  approximated  by 
u^(r)  according  to 

u  (r)  =  N  (r)  £  (4) 

where  £  is  the  vector  of  nodal  displacements  and  the  matrix  N  (r)  is  an 
"interpolation  operator." 

For  the  sake  of  illustrating  the  oftentimes  bewildering  finite  element 

notation,  it  is  convenient  to  use  the  simple  triangle  element  shown  in 

Figure  1.  Its  i^  node  is  at  (x^ ,  y^),  at  which  the  displacements  are 

(u  u  (i)).  Now  let 
x  y 

r  =  (x  y}H 


U  .  (ux  u/ 


iMuf'ufiufl  u<2>  i  u<3>  u<3>)H 

in  which  the  superscript  H  denotes  the  transpose. 

For  the  triangle  we  assume  the  linear  approximation 


ux  =  al  +  a2x  +  a3y 


Uy  =  a4  +  a  5X  +  a6y 


In  vector  notation 


u  =  D  a 


in  which  a  is  the  constant  vector 


a  =  (a,  a2  ^  a4 


and 
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J*\ _ 


'  1  x  y  o  o  o  " 
D  = 

.o  o  o  1  x  y  . 


It  is  elementary  to  derive  that 


where 


u  =  D  C  c 


C  = 


1  0  x-j  o  y,  o 

o  1  o  X]  o  yi 

1  o  *2  0  ^2  0 

o  1  o  x2  o  y2 

1  o  x3  o  y3  o 

o  1  o  x3  o  y3 


-1 


Returning  to  the  general  discussion,  the  true  strain  £  i 
be  written  as 

£  =  r  *u 

where  tT  is  a  kinematic  operator.  Applying  B'  to  £  furnishes 
approximation  as 

I  =  B'  *u 


=  B  C  £ 

where  B(r)  is  a  matrix. 

For  the  triangular  element  we  find 


1  =  {cxx  eyy  exy] 


from  which 


B  = 


0  1  0  0  0  0  1 

0  0  0  0  0  1 

Lo  0  h  1  1  l2 


The  true  stress  vector  £(r)  may  in  general  be  expressed 


n  an  element  may 


a  strain 


(5) 


as  a  functional 
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of  £,  with  the  actual  functional  form  determined  by  the  constitutive  model. 
Formally, 

a(r)  =  A(£,  r)  (6) 

The  approximate  stress  is  obtained  from 

o(r)  =  A(I,  r) 

=  A'  (_?,  r)  . 

In  the  triangle,  assuming  linear  isotropic  elasticity,  it  follows  that 

o  =  EBC  ' 

where  E  is  a  matrix  of  elastic  constants. 

For  equilibrium  of  an  element,  the  principle  of  virtual  work  may  be 
stated  in  terms  of  true  quantities  as 

J"  p  5udV  +  J'  £*  5zdV  -  J'  L  %dV  (7) 

V  V  s 

In  Equation  7,  V  is  the  element  volume  and  S  its  surface  area,  r  is  the 
traction  applied  to  the  element  boundary,  a  is  the  mass  density,  5( • )  is  the 
variational  operator,  and  the  superposed  dot  denotes  differentiation  with 
respect  to  time. 

We  assume  this  principle  also  applies  to  the  approximate  quantities: 

J'  p  6u_dV  +  J' iedV  =  J' <5udS  .  (8) 

V  V  s 

Hereafter,  the  overbars  designating  the  approximations  will  not  be  displayed. 
Upon  substituting  Equation  4,  the  right  hand  term  in  Equation  8  becomes 
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and  P  may  be  called  the  consistent  load  vector. 

For  the  inertial  term 

J'  ou^  5udV  = 

V 

-  *H  M  f.  - 

and  M  is  the  consistent  mass  matrix. 

The  second  term  on  the  left  hand  side  involves  the  constitutive  model. 
First  write  s  and  £  in  deviatoric  and  isotropic  components  as 

0  =  s  +  s  e 
£  =  e  +  e  _e 

where  £  is  the  vectorial  counterpart  of  the  Kronecker  tensor.  Elementary 
manipulation  leads  to 

0^  (5  £  —  sH  5  e  +  3  s  5  e  » 

From  the  constitutive  model 

£  =  2u  (e  -  ef) 
s  =  k  (e  -  ed) 

with  k  =  2v  +  3>..  Consequently, 

£^  o£  =  2u  eH  6e  +  3<  efe 

fH  H 

-  2m  e  Se  -  3k  e°6e 
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Equation  5  implies  that 


e  =  B'C 


,  Hr 
e  =  b  C 


(9a) 


(9b) 


where  the  matrix  B'  and  the  vector  b  are  easily  derived  when  B  is  specified. 
For  example,  in  the  triangle  element 

e  "  ^xx  +  £yy^3 


e  =  e  -  e 

XX  XX 


e  =  c  -  e 

yy  yy 


and 


ezz  =  -  e 


p  =  1  e  e  e  e  > 
—  xx  yy  zz  xy" 


'xy  ”xy 


Simple  algebra  leads  to 


b  =  {o  1/3  ooo  1/3} 


H 


B'  = 


o  2/3  o  o  o  -1/3 
o  -1/3  ooo  2/3 
o  -1/3  ooo  -1/3 
o  o  1/2  o  1/2  o 


Up  to  this  stage  in  our  development,  nothing  has  been  said  about  the 
distribution  of  and  in  an  element.  We  now  assume  that  they  are 
distributed  in  the  same  way  as  the  corresponding  parts  of  the  strain  tensor. 
Formally , 


e  =  B'C  3 
ed  =  bH  CY 


(10a) 

(10b) 


in  terms  of  new  unknown  vectors  3_  and  y_,  called  the  flow  and  damage  parameters. 
The  prime  in  Equation  10a  will  no  longer  be  displayed.  The  assumption 
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expressed  in  Equations  10a,  b  leads  us  to  call  the  present  relations  a 
consistent  inelastic  formulation. 

It  now  follows  that 


/oH  6  £  d  V  =  cVdt  +  CHK,£r 


h.  Jz.  d 

-  SKKf6£  -  YHKd5£ 


with 


f  ■  2  V  f 
T  V 

/ 


Kd  *  3t 


CHBHBCdV 


CHb*bHCdV 


and 


Ihe  matrix  Kg  given  by 


{ b*bn )  -j  j  =  b-jbj  . 


Ke  -  Kf  +  Kd 


is  nothing  but  the  ordinary  stiffness  matrix  of  linear  elasticity. 
The  equilibrium  relation  for  the  element  under  study  is  now: 


M  £  + 

(Kf  +  Kdk  =  P  +  Kf  S  + 

Kdl 

(ID 

In  the 

next  section 

we  use  the  constitutive 

model 

to  derive  equations 

governing  f 

and  y*  They 

will  have  the  general  form 

1  =  /  (l.  i,  kf ) 

(12a) 

C  f  f  Af 

kT  =  w  (£,  6,  kT) 

(12b) 

y  =  (t,  I.  kd) 

(12c) 

id  =  »d  (;,  V,  kd) 

( 1 2d ) 

where  z^,  w^,  z1^  and  w1^  are  material  functions.  More  concretely,  in  the  next 
section  the  constitutive  model  will  be  used  to  derive  the  material  functions 
in  Equation  12a-d. 
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B.  Finite  Element  Form  of  the  Constitutive  Equations 

The  constitutive  Equations  3a,  d,  e,  and  f  may  be  rewritten  as 


ef  •  £f  (s,  ef,  kf) 


kf  ■  hf  (s,  ef ,  kf) 


*  d  d  /  d  i  d  \ 
e  =  g  (s,  e  ,  k  ) 


kd  =  hd  (s,  ed,  kd) 


03a) 

(13b) 

(13c) 

(13d) 


Our  task  is  now  to  restate  the  constitutive  relations  in  terms  of  3_  and  y 
and  to  eliminate  dependence  on  r. 

From  the  constitutive  model  we  may  write 


e  =  BC  3 


(14a) 


s  =  2u  (e  -  ef)  -  2u  BC  (;  -  3) 


d  -  KHr 
e  =  b  C  y 


(14b) 

(14c) 


s  =  k  (e  -  ea) 


-  k  bMC  (;  -  v) 


(14d) 


Recall  that  B,  b^  and  C  may  depend  on  r. 

Note  that  k^  and  kd  appear  in  the  arguments  on  the  right  hand  side  of 
Equation  13a-d.  We  now  make  the  additional  assumption  that  they  may  be 

Af  A/j 

replaced  by  k  and  k  where  the  circumflexes  denote  the  element  volume 


averages: 


-  v  f  kf  dV 


(15a) 


kd  =  \  f  kd  dV 
v  V 


(15b) 
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From  Equations  13a-d,  14a-d  and  15a-b,  it  is  evident  that 
ef  =  £f  (2uBC(£  -  3),  BC8,  kf) 


=  (i>  3,  kf,  r) 


(16a) 


where  the  function  cf  is  defined  in  Equation  16a.  By  similar  argument, 

kf  =  pf  fq,  3,  kf,  r)  (16b) 

ed  =  qd  U,  x,  kd,  r)  (16c) 

kd  =  pd  (q,  x,  kd,  r)  (16d) 

We  seek  to  eliminate  dependence  on  £.  But  Equation  16a  implies  that 

f  CHBH  ef  dV  =  ( 2y ) ” 1 Kf  3 


=  / 


NHBH  £H  dV  . 


Therefore,  the  material  function  z_T  in  Equation  12a  is  given  by 

zf  (;,  s,  kf)  =  2uKf_1  f  CHNH  £f  dV 
-  -  -  T  V 


(17a) 


Similar  manipulations  serve  to  derive  w^,  zj ,  and  wd  in  Equation  12b-d: 


wf  -  V-1 


=  3k  K. 


/ 

V 

-1 


pf  dV 


/a 


qd  dV 


d  w-1  f  „d 
w  =  V  J  p 


dV 


(17b) 

(17c) 
( 1 7d) 


For  the  sake  of  illustration  we  consider  the  constitutive  relations 


e  -  nf  <  1  -  kf/Ff  >  (s.  -  cf  e  ) 


(18a) 
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Ff  =  f(i  -  cf  ef)H  (s  -  cf  ef)]1/2  (18b) 

ed  -  nd  <  1  -  kd/Fd  >  (s  -  cd  ed)  (18c) 


where  k^,  k^,  nd,  and  cd  are  material  constants.  These  relations 
were  previously  introduced  in  Reference  1. 

Applied  to  the  triangular  element,  the  relations  furnish 

S  =  vf  ff  (19a) 

in  which 

^  =  2  y  ?  -  (2p  +  cf)  S  (19b) 


vf  =  (2p  A)"1/2  [^H  Kf  4^]1/2 


(19c) 


where  A  is  the  area  of  the  element. 
For  damage. 


Y  = 


vd  'U 


with 


4^  =  <  -  (<  +  Cd)  Y 

vd  =  (3<A)"1/2  [^H  Kd  l^)VZ 


(20a) 

(20b) 

(20c) 


C.  Nodal  Continuity  and  Equilibrium 

The  previous  section  concerned  equilibrium  of  a  given  element.  Here  we 
consider  equilibrium  and  compatibility  of  the  assemblage  of  elements,  for 
example  the  triangles  shown  in  Figure  2.  For  this  purpose  it  is  adequate  to 
develop  force  balance  and  continuity  equations  holding  at  the  shared  nodes. 
Certain  modifications  of  the  single  element  relations  will  prove  convenient. 
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First,  some  additional  notation  is  needed.  The  quantities  _ 

■  and  now  refer  to  the  e"h  element,  for 

1  G  6  8 

example  the  second  element  in  Figure  2. 

Let  be  the  entry  of  j_'e^  referring  to  the  n"^1  node  and  the  mth 

( G ) 

direction.  For  instance,  since  ;  is  the  nodal  displacement  vector  for  the 


eth  element,  then  ' 

x 

^  is  the 

x-displ acement  of 

its  node 

having  n  , 

index.  In  reference  to  Figure  2, 

.(1)  =  ; ( 1  )  -  ( 1  ) 
—  ’X 

nun 

y 

(1)_(4) 

X 

(1 ) .  ( 4 ) 
y 

( 1 ) . (5) 

X 

H ) .  (5) 

V 

r ( 2)  =  -(2).(1) 
—  "X 

(2).(1 ) 
y 

(2). (2) 
"X 

(2), (2) 

'y 

(2) _ (4) 

X 

(2), (4)  ,H 

y 

,(3)  =  r(3)_(2) 

—  1  ’X 

(3), (2) 

y 

(3). (3) 

X 

(3). (3) 

y 

(3)  ,(4) 

'X 

( 3) . (4) tH 

y 

Continuity  of  displacements  implies  that 


OU1)  =  (2)r(l)  =  (1) 

’v  -v  uv 


(l).(l)  .  (2).0)  .  U(U 
y  y  y 


(2). (2)  ,  (3). (2)  =  (2) 

'’x  "x  x 


(2 )J2)  ,  (3). (2)  .  (1) 

y  y  y 


(3).  (3)  =  (3) 

'x  x 


(3) ,(3)  _  (3) 

y  y 


(21) 


(1),  (5)  ,  (5) 

"x  x 


0)r(3) 

y 


=  u 


(5) 


(Dr(4) 

”x 


(2). (4) 
“x 


(3). (4)  =  (4) 

"x  ux 


0)r(4) 

y 


(2). (4) 

y 


(3). (4)  =  (4) 

y  y 
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(  G  )  (  G  ) 

It  is  now  assumed  that  _  and  j_  '  are  continuous  in  the  same  sense 

{ G ) 

as  ^  .  This  is  not  implied  by  displacement  continuity,  but  it  is  expected 

to  assure  a  certain  degree  of  smoothness  in  the  distribution  of  the  flow  and 
damage  strains.  In  any  event,  the  alternative  would  appear  to  involve 
computing  a  prohibitive  number  of  inelastic  nodal  parameters. 


Suppose  for  instance  that  there  are  M  plane  triangular  elements  with  a 
total  of  N  nodes.  Under  the  present  assumption  there  are  2N  values  each  of 
and  to  compute.  But  otherwise  there  would  be  6M  values  of  and 

( G  ) 

to  determine,  and  M  is  nearly  twice  N.  Evidently,  the  inelastic 
smoothness  assumption  is  very  convenient  in  regard  to  computational  effort. 


Referring  to  Equation  21  and  Figure  2,  continuity  of  is  expressed 
as  follows: 


nun  ,  (2MD 


(1)?(1)  =  (2U1) 


(2)_(2)  m  (3h(2) 


( 2 )  ( 2 )  .  (3)  ,(2) 


0),(4)  .  (2)a(4)  ,  (3)  (4) 

^x  ~x  "X 


(22) 


(1 ) .(4)  .  (2)a(4)  .  (3), (4) 

y  y  y 


For 


Equation  22  holds  with  ganuia  substituted  everywhere  for  beta. 


Note,  however,  that  Equation  12a, c  must  now  be  modified  for  consistency 
with  the  nodal  continuity  of  and  v^e^.  Repeating  Equations  12a, c  in 
updated  notation 


£,e>  ■  4e>  (>>.  se,  4e>) 

•(e)  .  (e)  / .(e)  (e)  .(e), 

i  -  ±(j  u  *  i  »  K(j  i 


(23a) 


A  satisfactory  modification  is  to  replace  Equations  23a, b  with: 
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(e)  *  (n)  .  1_  y  (e)  (n)  -(e)  (e)  .  (eh 

-m  ‘  e_  4*  -  '  f  1 


n  e 


(e)-(n)  .  1  y  (e)  (n)  -(e)  (e)  .(e), 

o  *->  ZA  U  *  I  »  K  > 


en  7  "dm  i 


th 


where  en  is  the  number  of  elements  sharing  the  n  node. 

We  now  state  the  modified  relations  holding  at  the  fourth  node  in 

(41 

Figure  2.  First  define  ’  by 

.(4)  ,  1  ,(1L(4)  +  (2), (4)  +  (3)  (4), 

’’y  3  '  '’y  ’y  ‘’y 

and  by  virtue  of  displacement  continuity 

.(4)  .  (1  )„(4)  .  (2). (4)  _  (3)r(4) 

‘7  'y  "y  ”y 

(41  (41 

The  quantities  ‘  and  ‘  are  analogously  defined. 

Using  Equation  22  together  with  the  constitutive  model  represented 
Equations  12a-c,  it  follows  that 


;.(4)  . 


'4,  -!") 


(4) 


where 


*<4)  •  2,  44)  -  (2.  *  «f)  2<4> 

jj/e^  =  2p  -  (2y  +  cf ) 

v<«>  -!>*<'>  ]‘1/2  [(*<'))"  4°  *‘*)],/2 


i(e) 


th 


and  '  is  the  area  of  the  e  element 


For  damage  the  corresponding  relations  are 
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(24a) 


(24b) 


by 


(25a) 


4 
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Simply  stated,  the  sum  of  the  equivalent  reaction  forces  contributed  by  the 
e’ements  sharing  a  node  is  set  equal  to  the  external  force  aoplied  at  the 
node.  Equation  (26)  is  readily  extended  tc  more  general  situations. 

The  constitutive  relations,  for  examole  Equation  25  a,b,  together  with 
the  nodal  force  balance  relations  such  as  Equation  26  comprise  a  system  of 
ordinary  differential  equations  in  time.  Under  suitable  initial  conditions, 
the  eauations  may  be  integrated  numerically  to  furnish  the  nodal  displacement, 
flow  and  damage  parameters  as  functions  of  time. 

CONCLUSION 

The  finite  element  method  has  been  applied  to  a  constitutive  model 
describing  flow  and  damage  in  rapidly  loaded  structural  materials.  A  system 
of  ordinary  differential  equations  in  time  has  been  obtained  for  nodal 
displacement,  flow  and  damage  parameters.  The  formulation  is  "consistent" 
in  that  the  inelastic  strain  approximants  involve  the  same  interpolation 
operators  as  the  corresponding  parts  of  the  total  strain  approximant.  Certain 
interelement  continuity  conditions  are  imposed  on  the  flow  and  damage  strains. 
Numerical  results  will  be  reported  in  a  subsequent  article. 
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